Second-order shaped pulses for solid-state quantum computation 



OO 
O 

o 

(N 
G 



Oh' 



> 

O 

OO 

o 



X 



Leonid P. Pryadko 1 and Pinaki Sengupta 2 

1 Dept. of Physics, University of California, Riverside, CA 92521 
2 T-CNLS and MPA-NHMFL, Los Alamos National Laboratory, Los Alamos, NM 87545 

(Dated: June 16, 2008) 

We present the construction and detailed analysis of highly-optimized self-refocusing pulse shapes 
for several rotation angles. We characterize the constructed pulses by the coefficients appearing 
in the Magnus expansion up to second order. This allows a semi-analytical analysis of the perfor- 
mance of the constructed shapes in sequences and composite pulses by computing the corresponding 
leading-order error operators. Higher orders can be analyzed with the numerical technique suggested 
by us previously. We illustrate the technique by analyzing several composite pulses designed to pro- 
tect against pulse amplitude errors, and on decoupling sequences for potentially long chains of qubits 
with on-site and nearest-neighbor couplings. 

PACS numbers: PACS: 75.40.Gb, 75.40.Mg, 75.10.Jm, 75.30.Ds 



I. INTRODUCTION 

The implementation of .-quantum algorithms using 
NMR on molecules in liquidEJ, solidn, and liquid crystalstl 
has demonstrated in principle that pulse-based con- 
trol methods can be useful for quantum information 
processings (QIP). The technique has long been a sta- 
ple in NMR spectroscopy, where complex molecules like 
proteins are analyzed with the help of long sequences of 
precisely designed radiofrcquency pulsesO. Related tech- 
niques for coherent manipulation of many-body quantum 
systems have emerged as an important tool in many areas 
of science and technology. 

A useful quantum computer should contain hundreds, 
if not thousands of qubits. The only hope of scaling to 
such system sizes is with the help of multiple levels of 
quantum error correction (QEC). For this to work, the 
benefits due to each additional level of encoding should 
outweigh the corresponding overhead of additional errors. 
This leads to various threshold theoremsBQ, estimating 
the maximum error rate for which such concatenated 
encoding can be beneficial. The corresponding thresh- 
olds are rather stringent, meaning that for scalability one 
needs very accurate elementary gates. 

Even for relatively small n-body systems, the number 
of states scales exponentially with n, and the accuracy 
required for QIP is high. As demonstrated in several 
recent experiments in NMR QIP, required accuracy can 
be reached with the help of strongly-modulating pulses, 
where entire single- or few-qubit. gates are designed nu- 
merically for a given moleculcocl [also sec |I^|y],[l2[[13) . 
While the technique indeed offers unprecedentedly accu- 
rate, fast gates (which also helps to avoid relaxation), it 
obviously cannot be generalized to larger systems. 

In contrast, the traditional pulse and sequence design 
rely on the Magnus (cumulant) expansion!! The expan- 
sion is done around the evolution in .the ..applied con- 
trolling fields, while the chemical shiftstffijllj (resonant- 
frequency offsets) and inter-qubit couplings are treated 
perturbatively. The main advantage of the Magnus ex- 
pansion is its locality. Namely, when local qubit cou- 



plings are dominant, the control fields accurate to a given 
order can be designed by analyzing relatively small clus- 
ters. The results remain exact independent of the sys- 
tem size, or even in the limit of infinite system. One 
can thus characterize pulse-based method foij-designing 
control fields as scalable to large system sizes.E-3. 



A scheme to systematically construct high-order self- 
refocusing pulses and pulse sequences was developed by 
the authors4n-Rcf. [l^. Specifically, we constructed "soft" 
NMR-stylecMa second-order self-refocusing inversion (it) 
pulses and several high-order sequences based on such 
pulses for refocusing qubits arranged in spin chains with 
on-site chemical shifts and XXZ nearest-neighbor cou- 
pling. Theptnain technical advance which enabled the 
calculations^ was the efficient numerical algorithm for 
computing high-order terms of the Magnus expansion. 
The algorithm is based on the usual time-dependent per- 
turbation theory; the direct computation of multiple in- 
tegrals entering higher-order cumulants would be totally 
impossible. 



In this paper we present highly-optimized self- 
refocusing pulse shapes for rotation angles (f>o other than 
180°. For such pulses, we extend the results of Ref. 
and construct the analytical expansion of the evolution 
operator for an arbitrary coupled qubit. While the ex- 
pansion is more complicated than that for the inversion 
pulses with 4>q — 180°, to second order, it is still charac- 
terized by only three coefficients, two of which we sup- 
press by pulse shaping. This allows us to compute the 
error operators associated with a given control sequence 
semi-analytically, by evaluating the leading order terms 
in the corresponding products of the evolution operators. 
We illustrate the technique on several newly-constructed 
decoupling sequences for a chain of qubits with on-site 
and nearest-neighbor couplings, as well as with the com- 
posite pulses protecting against amplitude errors. 
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II. PROBLEM SETUP 
A. Dynamical Decoupling 

In principle, the simplest type of control pulses consist 
of short, intense bursts of coherent, resonant electromag- 
netic radiation, popularly known as "hard" or "bang- 
bang" pulses. In this limit, for the duration r of the 
pulse one can totally ignore all other couplings of the 
qubit(s). Then, a pulse sequence can be viewed as a 
series of free evolution intervals [unitary evolution opera- 
tors Uf(t) = exp(-itHs), where Hs is the system Hamil- 
tonian] intercalated with pulse operators. For example, 
with single-qubit control Hamiltonian, 



H, 



c 



(1) 



where cr^ 1 , /i ~ x,y,z arc the Pauli matrices for a-th 
qubit and VJf (t) are the corresponding control fields (or, 
more precisely, the envelopes of the control fields applied 
at the resonant frequency of the corresponding qubit), 
the pulse operator is the product of those for individual 

qubitS, P^Ua P a(4> { a \na), 

P a (0o , n) = cos — - in- a a any. (2) 

Here n is the unit vector that determines the spin rota- 
tion axis and 4>o is the corresponding angle, 



dtV»(t). 



The corresponding pulse algebra is straightforward. 
For example, for a single spin with the chemical shift 
Hamiltonian, 



H 8 = ±o>, 



(3) 



the sequence of two equally-spaced inversion pulses (</>o 
±7r) in x-direction is equivalent to a unitary, 



P(x,±)U f (t) P{-n,x)Uf{t) 
.tA 



(—ia x ) exp^— i^-a z ^j (icr x ) exp^— i— a' 
tA _\ / .tA 



/ tA 

ex P( +l ~5~ a ') ex P{-' l -^ a 



= i, 



(4) 



which is, of course, the fopnal identity behind the well- 
known spin-echo sequenccEil. 

In reality, the pulse duration r is always finite. Thus, 
during the pulse application, the rotation actually hap- 
pens around the axis determined by the sum of the 
system and control Hamiltonians; one gets finite-pulsc- 
duration errors. Generically, such errors scale linearly 
with pulse duration. These errors are especially danger- 
ous in systems where one cannot reduce the pulse du- 
ration indefinitely, e.g., because of the need for spectral 



addressing in homonuclcar NMR, or in order to avoid ex- 
citing -Levels outside the qubit state in Josephson phase 
qubitsc3. Yet, even in systems with optical qubit cou- 
pling where r could be in the sub-picosecond range, the 
finite-pulse-duration errors can be significant if one tar- 
gets the accuracy necessary for achieving the scalability 
thresholds^. 

The finite-pulse-duration errors can be significantly 
reduced by suppressing the leading-order error opera- 
tors. The errors would scale with a higher power of 
the pulse duration, which makes them much more man- 
ageable. This can be achieved either by designing spe- 
cial sequences (which typically doubles or quadruples the 
number of pulses), or by using specially designed self- 
refocusing pulse-shapesE-l Typically, the latter strategy 
is more efficientlia; besides, the self-refocusing pulses can 
often bejjp^d as drop-in replacements for corresponding 
(5-pulsc: 



B. Model 

We consider the following simplified Hamiltonian 

H{t) = H c (t)+H s + H v (t)+H <T , (5) 
with the first (main) term due to individual control fields, 
1 



Hc{t) 



2^ 



(6) 



where /i = x,y,z, are the usual Pauli matrices for 
the n-th qubit (spin) of the ID chain. The other terms 
include the native Hamiltonian of the system 



iJ« = - A fJ Tr fJ ' -I TP" n v n^ 1 



(7) 



describing the constant qubit couplings and the inter- 
actions between the qubits, and the coupling with the 
oscillator thermal bath, 



H a = \Y. n BZ<- (8) 



In Eq. (||), = A^{pi,qi) account for the possibility 
of a direct coupling of the controlling fields V£ with 
the bath variables q i} p i7 while P>£ = B£(pi,qi) describe 
the usual coupling of the spins with the oscillator bath. 
Already in the linear response approximation, the bath 
couplings (|^) produce a frequency-dependent renormal- 
ization of the control Hamiltonian Hc(t) [Eq. (0)], as 
well as the thermal bath heating via the dissipative part 
of the corresponding response function. Both effects be- 
come more of a problem with increased spectral width of 
the controlling signals V£ . In this work we do not spec- 
ify the explicit form of the coupling Hy(t). Instead, we 
minimize the spectral width of the constructed pulses. 
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While the Hamiltonian ([7]) is a generic spin Hamilto- 
nian, we will also consider specifically the Hamiltonian 
of XXZ model with additional on-site fields. 

H S XZ = \ E [ J n,n>«> +Jn,n> («'+«>)] 
(n,n') 

+^E A K- ( 9 ) 

C. Magnus expansion 

In a qubit-only system with the Hamiltonian 

H(t) = H c (t)+H s , (10) 

the effect of the applied fields is fully described by the 
evolution operator U(t), 

U(t) = Texp(-i J dt' #(*'))> (11) 

where T is the Dyson time ordering operator. 

For pulses with finite width, any desired unitary trans- 
formation can only be implemented approximately. A 
widely used framework to design pulses to effect a de- 
sired unitary transformation (or, equivalently, remove 
the effect nf-,undcsircd terms in the Hamiltonian) is 
the Maattuso expansion and the average Hamiltonian 
theoryEsEj. The expansion is done with respect to the 
evolution due to the control fields alone, 

[/b(t)= Texp(-*y ' dt'H c (t')), (12) 

by defining the system Hamiltonian in the interaction 
representation (the "rotating- frame Hamiltonian"), 

H s (t) = U^t)H s U (t). (13) 

For a periodic control field, H\{t + r c ) = H\{t), such 
that the zcroth-ordcr driven evolution is also periodic, 
Uo(t + t c ) = Uo(t), one has the following expansion in 
powers of r c : 

U(nr c ) = exp(-iHnT c ), (14) 
H = +BM +H (2 *> + ..., (15) 

where 

H {0) t c = I ' ' dtH u (16) 
Jo 

H {1) t c = -I f'dh f dt^H^Hx], (17) 
1 Jo Jo 

i rT e pt 3 pt 2 

H {2) r c = --J dt 3 J dt 2 j dti (18) 
x ([ff 3) [H^Ht]] + [H u [H2,H S ]]) ■ (19) 



Generally, the term If ( fc_1 ) contains a fc-fold integration 
of the commutators of the rotating-framc Hamiltonian 
Hi = Hs(U) at different time moments U and scales as 
||//^- ) ||t c oc j|T c i/g|| A; . Note that for small enough r c , the 
expansion parameter remains small even for long evolu- 
tion time. 

The Magnus expansion thus offers a basis for con- 
structing successful approximations towards the desired 
unitary evolution. With the simpler problem of decou- 
pling, the goal is to have no evolution. A K-th order 
refocusing sequence can be defined as that where there is 
no evolution to K-th order, that is, 

H (0) = = _ _ _ = ff(K-l) = Q (2Q) 

Respectively, at time nr c , the error in the unitary 
evolution operator would scale as \\U(nT c ) — 1| oc 
nllTcT^sll^" 1 " 1 , and the corresponding fidelity F(t) differs 
from unity by 

l-F{nr c )o,n 2 \\T c Hs\\ 2K+2 - (21) 

A crucial advantage of the cumulant expansion is 
that the cumulants do not contain the disconnected 
terms arising, ixpm different parts of the system (clus- 
ter theorcmO'Ea) . For an arbitrary lattice model of the 
form (^), with bonds representing the qubit interactions, 
the terms contributing to fc-th order can be represented 
graphically as connected clusters involving up to k lat- 
tice bonds; for a chain of qubit such clusters cannot have 
more than n = k + 1 vertices. Thus, to obtain the exact 
form of the expansion up to and including A'-th order, 
one needs to analyze all distinct chain clusters with up 
to K + 1 vertices. 

D. Time dependent perturbation theory 

While the Magnus expansion is conceptually straight- 
forward, it is cumbersome to implement and, most impor- 
tantly, the repeated integrations are very expensive com- 
putationally already at the second order, see Eq. jT^). 
An alternative strategy for evaluating high-order terms of 
the Magnus expansion was suggested by the present au- 
thors in Ref. [l8| Instead of working with the cumulants, 
the technique is based on the time-dependent perturba- 
tion theory (TDPT). The expansion is done around the 
non-perturbed evolution due to control fields alone, see 
Eq. ([l2"l). However, for actual computations, it is more 
convenient to use the differential equation 

U (t) = -iH c (t) U {t), U (0) = 1. (22) 

The slow evolution operator 

R(t)^U^ (t)U(t). (23) 

obeys the equation 

R(t) - -iH s (t)R(t), H s (t) = Ul{t) H s U (t), (24) 
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which can be iterated to construct the standard expan- 
sion R(t) = I + R 1 (t) + R 2 (t) + . . . in powers of (tH s ), 



R k (t) = -iH s (t)Rk-i(t), R (t) = I. 



(25) 



The successive terms can be evaluated by solving, at each 
step, a set of coupled first order ODE's simultaneously. 
For a finite system of n qubits and a given maximum or- 
der K of the expansion, one needs to solve Eqs. (|2^) and 
@ with 1 < k < K, These are total of (K + 1) coupled 
systems of first order ordinary differential equations for 
the 2 n x 2™ matrices Uq, R\, R 2 , ■ ■ ■ , Rk, and can be 
integrated efficiently. Computationally, this is a much 
less challenging job than that of evaluating repeated in- 
tegrals (|l7|) , (p"8|), and higher order terms. For a given 
system, solving the full set of equations (24) is simpler by 
a factor of at least (K + 1). However, it is the analysis of 
the perturbativc expansion that is the key for achieving 
the scalability of the results. 

Given the set of computed Rk(r c ) = Rk, the standard 
Magnus expansion can be readily obtained by taking the 
logarithm of the series, 



iH^T c = Ri, 



-iH^r r = R 2 - ^Rj, 



(26) 
(27) 



-iH^t c = R 3 -^[R 1 R 2 +R 2 R 1 ] +^Rl ...(28) 

Obviously, the order-if universal sclf-rcfocusing condi- 
tion (pO) is formally equivalent to 



Ri(t c ) = R 2 (t c ) 



= R k (t c ) = 0. 



(29) 



The matrices Rk in the latter condition arc much easier to 
evaluate numerically using Eqs. (|2^), (p5|). Importantly, 
the benefits of the cluster theorem are retained: to K-th 
order only clusters with up to K + 1 vertices need to be 
analyzed. 



where r is the pulse duration, fl = 2it/t is the corre- 
sponding angular frequency, and <fio is the requested ro- 
tation angle of the pulse. Note that the form d3f| ) guar- 
antees the symmetry of the pulse, V(t — t) = V(t). In 
addition, in order to reduce the spectral width of the 
control fields, we also constrained a certain number of 
derivatives of the function (|3f]) to vanish at t = and 
t = T, VW(0) = 0, Z = 0,1,...,2L-1. 

We implemented the computational algorithm de- 
scribed in the previous section using the standard fourth- 
order Runge-Kutta algorithm for solving coupled diL 
ferential equations (22), (|25|), and the GSL libraryEa 
for matrix operations. The coefficient optimization was 
done using a combination of simulated annealing and 
the steepest descent method. The target function for 
single-pulse optimization included the sum of the mag- 
nitudes squared of the matrix elements of the matrices 
Rk = Rk(r), k = I, . . . , K , as well as the weighed sum of 
the squares of the coefficients A m , 



K 



1/2 



f K = E tri? i^ 



(32) 



m— 1 



The second sum serves to provide some suppression of 
the higher Fourier harmonics of the pulse. In our simu- 
lations, the minimization was considered as having con- 
verged only after the first term evaluated to zero with nu- 
merical precision (typically, eight digits or better). For 
such a minimum to exist, the coefficient e in Eq. ( |32] ) 
should be sufficiently small (we used e = I0~ 4 ). 

For given pulse order K and the given number of addi- 
tional constraints L, there is a minimum number of har- 
monics M m i n (i<r, L) necessary for convergence. However, 
we found that the shapes obtained with M = M m - m (K, L) 
tend be over constrained and simply do not look nice. 
Our solution was to add one or two additional Fourier 
harmonics by increasing AI . 



III. PULSE DESIGN AND ANALYSIS 

A. Pulse design using TDPT 

The shapes of NMR-style one-dimensional pulsesQ, 
self-refocusing to a given order, can be found by ana- 
lyzing the single-spin dynamics with the system Hamil- 
tonian (^|) and the control Hamiltonian 

Hc(t) = \a x V{t). (30) 

Specifically, we encoded the trial pulse shapes in terms 
of their Fourier coefficients, 

V{t) = ^ +nJ2A m cos(mn(t-T/2)), (31) 



B. Pulse shapes 

Previously, in Ref. [T^, we gave the coefficients of the 
first-order self-refocusing (K = I) inversion (</>o = tt) 
pulse shapes Sl, as well as the second-order (K = 2) 
inversion shapes Ql, L = 1,2. Here L is the parameter 
that indicates the number of constraints at the ends of 
the interval: the value of the function and its derivatives 
up to (2L ~ l)st vanish at the ends of the interval [note 
that all odd derivatives are suppressed automatically due 
to the symmetry of the function, see Eq. j3f])] . 

In this work, we extend the list of constructed pulses 
to rotation angles </> = 1 0°, 20°, . . . , 180°. In Table @, 
we list the coefficients for the pulses used in the simula- 
tions. The coefficients for all of the constructed pulses 
are available upon request. 
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Ao 


Ax 


A 2 


A 3 


A A 


As 


Si(2tt) 


1.0 


-0.0237996956 


-0.6226198703 


-0.3535804341 






S 2 (2tv) 


1.0 


-0.0294359406 


-1.1741824154 


-0.2097531295 


0.4133714855 




Qi(2tt) 


1.0 


2.1406171699 


-2.3966480505 


-0.6474844418 


-0.0964846776 




Q 2 (2tt) 


1.0 


1.4818894659 


-2.6971749102 


-0.4384679067 


0.3434236044 


0.3103297466 


Si = Si{ir) 


0.5 


-1.2053193822 


0.4796467863 


0.2256725959 






S2 = S2 (?r) 


0.5 


-1.1950692860 


0.7841592117 


0.0737326786 


-0.1628226043 




Qi = Qi(tt) 


0.5 


-1.1374072085 


1.5774920785 


-0.6825355002 


-0.2575493698 




Q 2 S Q 2 (7T) 


0.5 


-1.0964843348 


1.5308987822 


-1.1472441408 


0.0025173181 


0.2103123753 


5i(tt/2) 


0.25 


-1.8963102551 


1.1337663752 


0.5125438801 






S 2 (tt/2) 


0.25 


-1.9049987341 


1.9858884053 


0.1063314501 


-0.4372211211 




Qi(tt/2) 


0.25 


-1.8948543589 


0.5873324062 


0.5970352560 


0.4604866969 




Q 2 W2) 


0.25 


-2.1145695246 


0.6415685732 


1.6854185871 


0.4511145740 


-0.9135322049 



TABLE I: Fourier coefficients for the constructed pulses, see Eq. (|3l|). Shapes Sl(4>o) and Ql(4>o) are the pulse shapes for 
rotation angle (f>o, respectively first (K — 1) and second (K — 2) order for the Hamiltonian (^). These shapes have 2L 
derivatives vanishing at the ends of the interval. 



C. Pulse shape analysis 



2. Leading-order average Hamiltonian 



The pulse shapes Ql{4>o), Sl{4>o) are constructed as 
first- or second-order sclf-rcfocusing pulses for the chem- 
ical shift Hamiltonian, Eq. (||). We would like, however, 
to have a universally good pulse shape that would work 
in most settings. To analyze the performance of the con- 
structed pulses in most general circumstances, we con- 
struct the Magnus expansion of the evolution operator 
for the most general system Hamiltonian, 



H S = A + A x a x 



A y ° v 



A z a z 



(33) 



where A p are the operators responsible for coupling with 
the outside worlds, and Aq is the external Hamiltonian. 
The analysis of the inversion pulses with (f>o = n appeared 
previously in Rcf. BO; here we extend it to (/>q 7^ 7T. 



1. Driven evolution 

The control Hamiltonian ( |30| ) alone [to zeroth order 
in Hs] produces the following unitary evolution operator 
[cf. Eq. (pi 



[7b (*) = e -^ (t)/2 , 4>{t)= / dt'V{t'). (34) 

Jo 

When acting on the spin operators, this is just a rotation, 
UQ(t)a v Uo(t) = <7 V cos 4>(t) — a z sin <fi(t). Consequently, 
the system Hamiltonian in the interaction representation 
has the form 

H s (t) = A a + a x A x + <r y (A y cos cj) + A z sin 0) 

+<7 z (A z cos(j)- A y su\(j)). (35) 



The zeroth order average Hamiltonian (|1 7| ) is just the 
average of Eq. ( p5| ) over the pulse duration. We assume 
V(t) represents a symmetric pulse, V(t p — t) = V x (t). 
Then, 4>(t) is antisymmetric, 4>(r p — t) = </> — <j)(t), where 
0o = 4>{ T p) is the overall notation angle. It is convenient 
to introduce the symmetrized rotation angle, 



<p(t) = 4>{t) - o /2, 



(36) 



such that <p(t p — t) = —ip(r p ). Then, the average of the 
sine over the pulse duration vanishes, (sin ip) = 0. This 
implies that the averages of the cosine and sine of the 
original rotation angle are 



c = (cos0) = cos(</>o/2)(cos</2), 
is = (sill 4>) = sin(0o/2) (cosy;), 



where 



If we denote 



(/(*)> = 



(cos ip) 



'p JO 



dtf(t). 



n dt 

/ — cos^(t), 
Jo T P 



(37) 
(38) 



(39) 



(40) 



then the zeroth order average Hamiltonian for a one- 
dimensional pulse becomes 



-ff (0) = A + a x A x -\ 
+<r z (^A z cos ■ 



<r y { A„ cos ■ 



A, sin 



2 J 



A y sin — 



(41) 



For the special case of the chemical-shift Hamiltonian 



we have Ao 
gives 



A x 



= y — I (j~ cos 

2 2 
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V,(t)/(2)l) 
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gale error: 2. 17743e-10 './' •-...•* 
0.2 0.4 0.6 




(c) 

FIG. 1: A single-spin second-order inversion (tt) pulse Qi(n). 
(a) Pulse profile over a complete period, (b) Evolution of the 
spin beginning with s z = 1. (c) The power spectrum of the 
pulse. The vertical lines denote the location of the harmonics. 
As seen, the spectral weight is almost entirely confined to 
u> < 5uiq. 



(b) 



gale error: 9.22685e-10 '.„•' 

0.2 0.4 0.6 




(c) 

FIG. 2: As in Fig. [j], but for the single-spin second-order it/2 
pulse Qi(tt/2). 



Clearly, the lst-order self-rcfocusing condition corre- 
sponds to u — 0. For such pulses the full zeroth-order 
average Hamiltonian is given just by the two first terms 
in Eq. ©. 



3. lst-order average Hamiltonian 

The lst-order average Hamiltonian (|l7]) is given by a 
double integral of the commutator of the system Hamil- 
tonian in the interaction representation (|3tj) evaluated 
at two different times. We note that every term in 
Eq. (|35| ) can be classified as either time-independent [pro- 
portional to e(t) = 1], proportional to c{t) = cos(f>(t), or 
to s(t) = s'm(f>(t). Therefore, most generally, the second- 
order terms in the evolution operator can contain the 




FIG. 3: (color online) Pulse shapes for (f>o = n/2. Solid 
lines represent Ql (7r/2), dashed lines correspond to Sl(tt/2). 
Pulse shapes with L = 1 are drawn with thin blue lines, while 
those with L = 2 are drawn with thick red lines. The black 
doted line shows the Gaussian shape Goio(t/2). 
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FIG. 4: (color online) As in Fig. g but for the inversion pulses, 
4>o — n. Note that the 1st order pulses (Sl(tt), dashed lines) 
actually have a smaller power than the Gaussian pulse. 




FIG. 5: (color online) As in Fig. y but for the pulses with 
4>o — 2n. The pulse shapes appear to resemble those of a pair 
of consecutive n pulses. Second-order pulses happen to have 
the smallest power. 

following nine integrals, 

ee = (l'l) = E (l'cosc/)), es = (l'sin^), 
ce = (cos </>' 1) , se = (sin (/>' 1) , 
I 



cc = (cos 4>' cos 4>) , cs = (cos cj>' sin (j>) , 

sc = (sin 4>' cos 0) , ss = (sin <p' sin 0) , (43) 

where we used the notation 

</(0')<?(0)> = -4 C dt'f{m) f dtg(<f>{t))> (44) 
r Jo Jo 

and 1 = e(i) or V = e(i') indicate an identity factor at 
the corresponding position of the average. However, be- 
cause of the commutator structure in Eq. (|T7|), only the 
following antisymmetric combinations appear in the ex- 
pression for the corresponding term in the average Hamil- 
tonian theory, H^ 1 ', 



a sc — cs ec — ce es — se 

2=^—> Cc = ^2~ ' Cs = ^^' (45) 

where 

a = ^ [ V dt 1 [ dt sm(<j>{t') - <j){t)) (46) 
2 r p Jo Jo 

In fact, the coefficients Cc and Cs can be reduced further, 

J- J- • 00 j- ,00 / a r-r\ 

Cc = Csm— , Cs = -Ccosy, (47) 

where [see Eq. ( ^6|) for the definition of <p(i)] 
f T " dt' f t 1\ . , , 

Thus, to second order, the average Hamiltonian of a sym- 
metric angle-0o one-dimensional pulse is determined by 
only three dimcnsionless coefficients, v, a, and £, see 
Eqs. (^5|) , (HI), and (|§|). These coefficients contain all 
the relevant information about the shape of the pulse. 

An explicit calculation of the lst-order average Hamil- 
tonian gives 



ffW = aT p (i[A z ,A y }-a x (Al + A 2 z )) 

+ Ct p cos y (a v (i [A z , Ao] +{ A x , A y } ) + o z (i [A , A y ] + { A x .A z } )) 

- ^ sin ^-(a«(i{A v , Ao] - {A X ,A Z }) + o z (i[A z , A ] + {A x , A y })) . 



For the Hamiltonian (^), the terms with C disappear, and 
we have, simply 

H U=-aa x ^. (49) 

Thus, the second-order self-rcfocusing pulses have both 
v = and a = 0. 



The actual parameters for the pulses with <pQ = 7r/2, 
7r, and 2tt are listed in TAB. [fi|. 
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pulse 


00 


w = (cos if) 


a 


c 


(f> S(t - t p /2) 


4>o 


00 
COS 

2 


sin <f>o 
4 


1 00 

- sin — 

4 2 


2*(*-t,/2) 


tt/2 


^2/2 


1/4 


V2/8 


G . 05 [90] 


tt/2 


0.730111 


0.398519 


0.175999 


G .i(90) 


tt/2 


0.753116 


0.420275 


0.173665 


Si (90) 


n/2 





-0.013067 


0.198719 


02(90) 


7T/2 





—0.0294665 


0.182109 


Ql(90) 


7f/2 








0.202067 


Q 2 (90) 


7t/2 








0.161658 


7r(t - r p /2) 


7T 








1/4 


Go. 05 (180) 


7T 


0.0744894 


0.0377451 


0.249476 


G .i(180) 


7T 


0.148979 


0.0764911 


0.247905 


Si (180) 


7T 





0.0332661 


0.238227 


S2(180) 


7T 





0.0250318 


0.241378 


Qi(180) 


7T 








0.239888 


Q 2 (180) 


7T 








0.242209 


2-ir(t - t p /2) 


2tt 


-1 








Go. 05 (360) 


2tt 


-0.896959 


0.402852 


0.00291436 


G .i(360) 


2tt 


-0.793918 


0.317488 


0.0116577 


Si (360) 


2tt 





0.0739621 


0.113233 


S 2 (360) 


2tt 





0.0612747 


0.0811486 


Qi(360) 


2tt 








0.00403872 


g 2 (360) 


2tt 








0.00734526 



TABLE II: Parameters of several common pulse shapes. The 
first line represents the "hard" ^-function pulse applied at the 
center of the interval of duration r, Gooi denotes the Gaussian 
pulse with the width 0.01r p , while S„ and Q n denote the 1st 
and 2nd-order self-refocusing pulses from Tab. @. 



IV. OPEN SYSTEMS 

In this work we concentrate on the performance of 
high-order pulses and pulse sequences in closed quantum 
systems. However, it turns out that such sequences also 
remain efficient in op(»,|Sjfstems, in the presence of low- 
frequency bath modes£SE3. 

The analysis is done in general form with the help of 
an assumption that the bath couplings have the same 
form as the existing terms in the system Hamiltonian (Q) , 
which are assumed to be suppressed to order K = 1 
or K = 2. The bath modes are assumed to be low- 
frequency; in addition to the expansion in powers of the 
corresponding couplings, one needs a low frequency ex- 
pansion in powers of the adiabaticity parameter t c /tq, 
where r c is the decoupling cycle duration and To is the 
bath correlation time. 

With K = 1 decoupling, the effect in the open sys- 
tem is a suppression of direct decay (Ti) processes, as 
well as the reduction of the dephasing rate (T2) by the 
factor of order of the adiabaticity parameter t c /tq. The 
former result can be understood by analyzing the spec- 



tral properties of the driven systemEilE3, while the latter 
can be viewed as due to a reduction of the time step for 
phase diffusion. With second-order decoupling, K = 2, 
the decoherence rate is additionally suppressed, and with 
time-reversal invariant bath coupling all orders of the ex- 
pansion in powers of adiabaticity parameter may vanish, 
in which case the leading-order dephasing term becomes 
exponentially small and dephasing would likely be deter- 
mined by terms of higher order in bath coupling. Along 
with the decoherence rates characterizing the exponen- 
tial decay of quantum correlations with time, the corre- 
sponding prefactor, which determines the ".visibility" (or 
"initial decoherence"E3) , was also analyzed^. While for 
generic refocusing sequences with K > 1 the initial de- 
coherence is quadratic in t c and does not scale with the 
thermal bath correlation time To, for symmetric pulse se- 
quences it is reduced by an additional power of the adia- 
baticity parameter (t/tq). These results were originally 
derived for a generic featureless bath, but they also hold 
in a vicinity of a sharp resonance as long as the effective 
(i.e., renormalized as in the average Hamiltonian) cou- 
pling jta the corresponding mode is small compared to its 
widthE3. 



V. APPLICATION EXAMPLES. 

A. Decoupling sequences for a chain of qubits 

Decoupling sequences are designed to prevent quantum 
evolution from happening. Thus, we want to construct a 
sequence such that the resulting evolution operator over 
the period t c is identity, U(t c ) = 1. We illustrate the 
scalability of dynamical decoupling to large system sizes 
by considering linear chains of qubits with cither Ising or 
XXZ n.n. random-valued couplings [only J* n+1 or both 

J^n+l and Jln+1 = Jn,n+1 in Ec h ©L P luS thc loCal 

fields cither along z axis or in arbitrary direction [A^ ^ 
or 7^ for fj, = x,y,z in Eq. <uu]. 

With such a system Hamiltonian, zcroth-order average 
Hamiltonian ([l6]) contains only individual qubits or pairs 
of neighboring qubits, the largest clusters contributing 
to the lst-order average Hamiltonian (O) originate from 
two bonds sharing a site (three qubits]7 and in general 
contains terms spanning contiguous clusters of up 
to n + 1 bonds, that is, n + 2 qubits. Thus, to design 
a K-th order decoupling sequence, one needs to consider 
individual clusters of up to K + 1 qubits. 

With nearest-neighbor and local couplings only, the de- 
coupling can be implemented by simultaneously applying 
pulses on either odd or even sublattice. We note that in 
our setup there is no gap between subsequent pulses, the 
pulses follow back to back with the repetition period r. 
The system is "focused" at the end of each cycle consist- 
ing of several pulses of length r. Such a scheme with a 
common "clock" time r is convenient, e.g., for parallel 
execution of quantum gates in different parts of the sys- 
tem. For each qubit, various pulses (or intervals of no 
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signal) can be executed in sequence. 

In this work we consider the following two sequences 
from Ref. [l8[ 4 = X 1 Y 2 X 1 Y 2 and its symmetrized ver- 
sion 8 = X1Y2X1Y2Y2X1Y2X1, which provide universal 
refocusing of the couplings between the sublattices, and 
also suppress the on-site chemical shifts A*. Here, Xi 
is a 7r x pulse simultaneously applied on all odd sites, Y2 



tt- v pulse applied on all even sites, etc. 



is a (-7r)i 

These sequences are "best" sequences at given length for 
all pulse shapes found by exhaustive search (high-order 
sequencesoEil equivalent for hard pulses do not necessar- 
ily have equal orders here). The fact that such a brute- 
force optimization approach works is entirely due to the 
efficiency of the numerical method. 

In addition, we constructed two longer sequences, 
16 = X1Y2Y1 OX1X2 Y1OX1Y2Y1OX1X2Y1O, and its sym- 
metrized version 32, constructed by running the sequence 
16 first directly and then in reverse order. Here denotes 
zero pulse, an empty interval of duration r. These two 
sequences provide universal decoupling both for any cou- 
plings between the sublattices and for arbitrary on-site 
fields (A£ £ 0). 

In addition to the system Hamiltonian, the effective- 
ness of a sequence application depends on the quality of 
the pulses. In Table III, we list orders of the sequences 
when applied with different pulse shapes, computed us- 
ing the numerical time -dependent perturbation theory as 
described in sec. 



II D 



The term i?/t(r c ) was considered 
to be zero if its norm vanished with numerical precision, 
typically or better, compared to typical values of 

order one for orders where Rk{r c ) 7^ 0. The orders K do 
not depend on the chain length; wc verified this state- 
ment on chains up to n = 7 qubits. Also, the computed 
orders are the same for all self-refocusing pulse shapes of 
particular order; we believe that the results will remain 
valid for other symmetric pulse shapes of the same order 
as indicated in the 1st column of Table III. 



B. Error scaling 

We illustrate the predicted power laws in Fig. ^, where 
the average infidelity (A4) is computed for different ra- 
tios of t/r, where t is the fixed evolution time and the 
pulse duration r was reduced to accommodate a differ- 
ent number of decoupling cycles. The simulation is done 
for chains of n = 4 qubits with randomly chosen but 
fixed parameters corresponding to different chain mod- 
els as indicated. The steepest lines correspond to largest 
order K of the sequence decoupling order. For symmet- 
ric sequence 8 with Ising chain, K — 2 for Gaussian 
pulses, Fig. ^(a), K = 4 for lst-order pulses, Fig. ||(b), 
and K — 6 for 2nd-order pulses, Fig. [|(c). The cor- 
responding infidelities for fixed evolution time scale as 
{Jzt) 4 , oc (J z t) 8 , and oc (J,t) 12 . Larger values of 
K can improve accuracy by orders of magnitude, or, at 
fixed required fidelity, substantially reduce the number 
of decoupling cycles. 



pulse 


model 
sequence 


Ising 


I+A? 


xxz 


XXZ+Af 


XXZ+Ai 


Ql, 


4 


5 


2 


1 


1 





all K = 2 


8 


6 


3 


2 


2 





pulses 


16 


2 


2 


1 


1 


1 


(w = a = 0) 


32 


3 


3 


2 


2 


2 


Sl, Herm jwj 


4 


3 


1 


1 


1 





all K = 1 


8 


4 


1 


1 


1 





pulses 


16 


1 


1 


1 


1 


1 


= 0) 


32 


1 


1 


1 


1 


1 


!^ 

Gauss Q 


4 


1 














8 


2 


1 


1 


1 







16 



















32 


1 


1 


1 


1 


1 



TABLE III: Order K for several decoupling sequences used 
with different pulse shapes upon different spin chains with 
nearest-neighbor and local couplings. Order K means that 
the first non-zero term in the average Hamiltonian ([wj) is 
H^ K \ so that for small enough r the mismatch in the unitary 
evolution operator after n decoupling cycles (evolution time 
t — nr c ) scales as \\U — 1|| oc tr^ , and the corresponding 
infidelity 1 — F oc t 2 r^ lc . Sequence 8 a sequence of 8 pulses 
applied intermittently on odd or even sublattices, see text for 
actual definitions. 



We saw that with order-A decoupling in multi-qubit 
systems with local couplings, the decoupling error oper- 
ators can be represented as connected clusters of up to 
K + 1 bonds. For a linear chain, these involve up to 
K + 2 qubits, and the number of such operators scales 
linearly with the total number n of qubits, as long as 
n > K + 2. In an n-qubit system, each of such operators 
can be written as an outer product of the cluster con- 
tribution, and the identity operators for the remaining 
qubits. As a result, the square of the Frobenius norm 
of the error operators scales linearly with the size of the 
Hilbert space, that is, exponentially with the number of 
qubits. However, this exponential scaling is suppressed 
when we compute the infidelity [see Eq. (A4)], so that the 
infidelity scales only linearly with the number of clus- 
ters, that is, linearly with the number of qubits. The 
same scaling with the system size is expected in higher 
dimensional arrangements of qubits (planar, 3D). 

We illustrate the scaling of decoupling errors with the 
qubit number n in Fig. [?]. The plots show the scaling of 
the average infidelity at the end of the interval in Fig. ^ 
and other data with the chain length n. 



C. Composite pulses 

Composite pulses are, in fact, pulse sequences designed 
to replace a single pulse and specially designed to com- 
pensate for some particular systematic errors, includ- 
ing off-rcson.ance. flPFlUf&tip 1 ]!! If| u l sc amplitude, and pulse 
phase errorsMOGa1^ oElc3cjCJ. 
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10"' 



XXZ+A 2 8K=1 □ 

lsing+A z 32 K=1 * 
lsing+A z 8 K=1 

•Jsing 8 K=2 + 



(a) 



(b) 
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10 

Jiff 
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10 
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S, n=4 XXZ+A^ 1 32 K=1 
XXZ+A Z 8K=1 
lsing+A z 32 K=1 
lsing+A z 8 K=1 
Tsing 8 K=4 
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Q, n=4 XXZ+A^ 1 32 K=2 
XXZ+A Z 8K=2 
lsing+A z 32 K=3 
lsing+A z 8 K=3 
Ising 8 K=6 



(c) 



FIG. 6: (color online) Illustration of decoupling accuracy with 
sequences 8 and 32 for chains of n = 4 qubits with different 
couplings as indic ated on the plots. The plots show average 
infidelity [see Eq. ( A4 )] computed at fixed time t as the pulse 
duration r was reduced to accommodate a different number 
of sequences. The values of model parameters were randomly 
chosen and remained the same for all simulations. Symbols 
are the data points, li nes are the single-parameter fits of the 
mismatch 8 [see Eq. (A3)] to 8 — br K , where the values of 



K indicated on the plots correspond to those in Tab. III. (a) 
Gaussian pulses; (b) lst-order pulses Si; (c) 2nd-order pulses 



The off-resonance errors appear when the carrier fre- 
quency of the applied pulse is off the transition frequency 
between the |0) and |1) state of a qubit. In the rotating 
reference frame this is equivalent to a non-zero chemical 
shift Hamiltonian (^|) . with A equal to the corresponding 
frequency bias. We note that our 1st and 2nd-order self- 
refocusing pulse shapes already offer a degree of stability 
against such errors. 

For this reason we concentrate on the pulse amplitude 
errors, where the correct pulse shape is applied with the 




FIG. 7: (color online) Scaling of the infidelity 1 — F at 
t/r = 128 with the chain length n for a particular realiza- 
tion of an XXZ chain with on-site disorder Af, decoupling 
sequence 8, pulse shapes as indicated. (Data for the pulse 
Goio divided by 10 to fit with the other data.) While the uni- 
tary matrix mismatch 8 2 [see Eq. (A3)] grows exponentially 



with the chain length n, 5 oc 2 n , the leading-order contri- 
bution to the corresponding infidelity (1 — F) represents the 
probability of error in one of the clusters, and it scales only 
linearly, as also seen in the plots. 



wrong amplitude, producing an incorrect rotation angle 
0o 7^ 4>o- Note that no one-dimensional pulse shaping 
can compensate for this kind of errors, since the mod- 
ified rotation angle is simply proportional to the pulse 
amplitude, 0o = (1 + f)4>a- 

On the other hand, one can expect that the pulse am- 
plitude offset / remains the same for all the pulses applied 
at a particular frequency. This uniformity is utilized in 
several composite pulses designed so that the net rotation 
would be insensitive to such uniform errors. 



1. SCROFULOUS 

The three-pulse sequence SCROFULOUSLiLis, based on 
the sequence originally proposed by TyckcOCJ. Particu- 
larly, an improved 7r pulse is obtained by applying three ti 
pulses, at 60°, 300°, and again at 60°, or just ttqq ^300^60- 
In the case of ideal <5-pulses, the resulting pulse compen- 
sates for pulse amplitude errors to linear order. With 
finite-width shaped pulses, an additional error is gener- 
ated due to the presence of the system Hamiltonian. In 
particular, for the chemical-shift Hamiltonian (||), the ex- 
pansion of a unitary operator applied along x axis has the 
form [see Eqs. @), @] 



U x = cos — - ia x sin — 

-1— y~ (cos</> er z - G y sin 0o ) 
, t 2 AV (. . 0o 00 

H ^ l&x Sm -T. C0 S -TT 

8 V 2 2 

r 2 A 2 a ( . 9 
H ^ — I icr x cos — + sm — 



^+ S in^ )+0(ry. (50) 
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Combining the corresponding expressions appropriately 
rotated in the x-y plane, and expanding the result to 
quadratic power in the relative amplitude offset /, we 
obtain for the composite pulse 7T6o ^300^60, 



CR 



-ia x + ir Aver, 



« Z CTi, 



(51) 



where v = v(f) = v + v'f + C(/ 2 ) is the parameter v 



Eq. (40)] but for the pulse with rescaled amplitude, and 



the further terms are of order fvrA, v r A , afr A . 
Clearly, with Gaussian or other pulse shape such that 
v 7^ 0, the error is linear in tA and quadratic in the 
amplitude shift / (although generally there will also be 
a cross-term cx /rA). This situation is illustrated in 
Fig. ||(a), where the average infidelity is plotted for the 
SCROFULOUS sequence with pulses G010 on the plane 
At-/. The region for 1 — F = 1CP 5 is a narrow vertical 
line which corresponds to great sensitivity to frequency 
shift. With lst-order self-refocusing pulses such that v = 
0, v(f) cx /, the error is dominated by the term cx /tA. 
The corresponding region corresponds to a diamond-like 
shape in the center of Fig. ||(b). We have also generated 
self-refocusing pulse shapes such that both v = and 
u' = 0. Then, by symmetry, v" = 0, and generically 
v cx / 3 . Then, for lst-order pulses, a^O, the errors are 
dominated by the term omitted in Eq. (fill); they scale 
as C(/ 2 ), 0(v"'Af 3 ), 0(oA 2 ), while for second-order 
pulses the last two terms become 0(A 3 ) 0(a' f 2 A 2 ). 
Plots for such shapes are shown in Figs. g(c) and ||(d) 
respectively; the result of improved pulse stability is a 
much wider region of high fidelity. 



2. BB\ and related pulses 

A longer but more accurate composite pulse known as 
BBi was originally proposed by WimperistlH. For tar- 
get angle 6 = tt, the pulse can be written as BB' 1 ^ = 
7ro7i>(27r) 3 07r0, where 4> = — cos _1 (— 1/4) sa 104.5°. For 
ideal <5-pulses, this cancels errors of both 1st and 2nd or- 
der in the relative pulse amplitude bias /. A related sym- 



metrized sequence BB 
was proposed in Ref. 38 



(CLJ) 



(7r/2) O 7r0(27r) 3 07r0(7r/2) o 
[see also Ref. £59|; because of 
the symmetry it leads to some additional error cancella- 
tion at higher order. With shaped pulses, we have also 
analyzed variants of these sequences with the 2tt pulses 
replaced by two tt pulses, BB^ ' = TroTr^^^^^s^TTcj, 

and Bb[ CLJ ) = (71/2)0^(^)30(^)30^(^/2)0. 

Computing the products of versions of Eq. (^0|) ap- 
propriately rotated in the x-y plane, for on-resonancc 
application of any version of the BBi sequence, 



17, 



(A=0) 
BBi 



L^ { ^ ll5 ^ az) + (D{n (52) 
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FIG. 8: (color online) Contour plots of the average fidelity 
for the composite pulse SCROFULOUS 7^0^300^60 with (a) 
Gaussian pulses G010; (b) 1-st order self-refocusing pulses Si 
(the plots for pulses Qi look similarly but symmetric with 
respect to horizontal axis); (c) lst-order pulse with ampli- 
tude correction v = v' — v" = 0; (d) 2nd-order pulse with 
amplitude correction, v — v' = v" — a = 0. The axes are 
the relative frequency mismatch tA and the relative pulse 
amplitude (1 + /), see text. 



We note that to achieve the level of infidelity of, say, 1 — 
F = 10~ 4 , the frequency mismatch should satisfy |/| < 
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0.136, compared with |/| < 0.090 for the SCROFULOUS 
sequence. Thus, even though order of the BBi family of 
composite pulses is higher (and thus, for small |/| their 
performance is much better asymptotically), at this level 
of infidelity their performance is comparable. 

We now turn to off-resonance correction terms which 
differ between implementations of the BBi sequence. In 
particular, with generic pulses such that v ^ 0, already 
at / = 0, all of these sequences acquire linear correc- 
tions scaling with rA. For example, the expansion of the 
sequence BB 1 ^ at / = can be written as 



Tj(W,f=o) _ _ l(T 

U BB 1 — la i 

T 2 A 2 



+ T (2icr x v 2 + 6v 2 t-'i + a + iVlEcr z (a - 2viv 2 )^ 

+0(A 3 r 3 ), (53) 

where a = 2a\ — a 2 , v 2 = vf+v 2 , and the parameters v±, 
ot\ and v 2 , a 2 correspond to the 7r and 2tt pulses respec- 
tively. Notice that the second-order coefficients oli enter 
only in the combination 2ai — a 2 . Not surprisingly, if we 
replace the 2-7r-pulsc with two 7r pulses, the coefficients a 
cancel out, 



(W',f=0) _ 



BBi 



tAvi 

+ l T^Oz 



+ i- 



(54) 



thus second-order accuracy can be obtained already with 
1st order pulses, vi = 0. Now, when both the ampli- 
tude and the resonant frequency bias are present, / =/= 
and A ^ 0, there is an additional source of error due 
to dependence of the pulse parameters on the ampli- 
tude, v, a — ► £>(/), <5(/). For regular self-refocusing 
pulses, v(f) = v'f + 0(f 2 ), and the 1st order terms in 
Eqs. (53), ( p4[ ) and their analogs for the other variants of 
the BBi pulse dominate the error cx /rA. Respectively, 
the average infidelity scales as oc (/rA) 2 , resulting in 
characteristic diamond-like shape on the contour plots of 
infidelity, see Figs. |9](a,b), |l0|(a). 

With specially designed pulses such that both v = 
and v' = 0, due to pulse symmetry also v" = 0, so 
that v(f) = v"'f 3 + 0(f i ). Then, for sequences other 

than BB^ \ with lst-order pulses the error is domi- 
nated by the terms quadratic in rA due to coefficients 
a, [cf. Eq. ( |53|) . As a result, the high-fidelity regions 
in Figs. p|(c), p-O(b), and [ll](a) have much more rounded 
shape. With 2nd-order pulses with amplitude correction, 
such that on = but a[ ^ 0, the leading-order error 
scales as C(/r 2 A 2 ), which extends the high-fidelity re- 
gions out to larger values of r A in a characteristic "smile" 

pattern. For the sequence BB^ W ^ these terms cancel out, 
and the leading-order error term comes from the non-zero 
v'", the errors scale as oc / 3 rA, 



0(r 3 A 3 ). (55) 



in a z . The result is a somewhat skewed in the center 
high-fidelity region widely stretched horizontally. 



D. Stability of decoupling against amplitude errors 

We now return to the problem of decoupling for a chain 
of qubits, but now consider the effect of the amplitude 
errors. This needs a separate study since single-qubit 
errors with composite pulses have structure different from 
those due to, say, finite pulse width. 

In Fig. |lj(a) we present the results of simulations 
for a particular 4-qubit Ising chain with on-site chemi- 
cal potential shifts Af over time interval t = 128r (ex- 
actly the same parameters as in Figs. [| 0). Specifi- 
cally, we plot the infidelity (1 — F) in units of 10~ 4 , 
as a function of the fractional pulse amplitude 1 + 
/. In the first half of the symmetric sequence 8 = 
XiY 2 XiY2 F 2 XiY 2 Xi, we used the original BB^° pulse 
X — * ^0^(2^)30^0 for Tlx, rotated appropriately to im- 
plement Y -> 7r 7r / 2 7r0 +7r / 2 (27r) 3 +T /27r0 +7r /2 as well as 
X = Xjr, while in the second part of the sequence 
we used the same decomposition but backwards, e.g., 
X -> 7r (27r)307r 7r o . Here <j> = -cos _1 (-l/4) w 104.5°. 
With no amplitude mismatch, thus constructed decou- 
pling sequence has order K = 1 with lst-order pulses, 
and an almost-2nd order "1*" with 2nd-order pulses, 
with the norm of the 2nd-order term reduced by three 
orders of magnitude compared to lst-order pulses. These 
should be compared with K = 1 and K = 3 respectively 



This error has the same order in / as that in Eq. (52), 
and it can compensate or increase the contribution linear 



for the regular 8-pulse sequence [Tab. III]. 

It is seen from Fig. |l2|(a) that regular 1st- and 2nd- 
order pulses Si and Qi perform in the BBi-corrected 
composite sequence almost as well as the self-refocusing 
pulses with additional amplitude protection (lst-order Li 
and 2nd-order L3). In addition, the 2nd-order pulses 
with amplitude protection (L3) work almost as well in the 
regular 8-pulsc sequence. All of these allow to achieve the 
infidelity level of 10 -4 at amplitude mismatch of 3% or 
higher (up to 6% with pulses L3 in the sequence 8BB1). 
With the regular 8-pulse sequence, the region of high 
fidelity shrinks substantially for pulses Qi and Li, and it 
all but disappears for the pulse Si [the coefficient a for 
the pulse Li happens to be a few times smaller than that 
for Si}. 

The results of analogous calculation for a particular 
4-qubit XXZ chain with on-site chemical potential shifts 
Af over time interval t = 128r (exactly the same parame- 
ters as in Figs. ||, [?]) are presented in Fig. |lj(b). However, 
it turns out that the BBi composite pulses lose accuracy 
when used with XXZ chain. Even in the absence of am- 
plitude errors, there are linear errors in tJ 1 - [the zeroth 
order average Hamiltonian is non-zero]. Thus, we only 
present the results for the regular 8-pulse sequence. With 
this sequence, the decoupling orders for 1st- and 2nd- 
order pulses are K — 1 and K = 2 respectively [Tab. III]. 



Compared with Ising-only couplings, with the particular 
parameters chosen, this increases the infidelity by some 
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FIG. 9: (color online) Contour plots of the average infidelity 
1 — F for the composite pulse BB^ W ' [710^^(271)3^71,1,] with (a) 
lst-order self-refocusing pulses Si; (b) 2nd-order pulses Qi; 
(c) lst-order pulses with amplitude correction; (d) 2nd-order 
pulses with amplitude correction. The axes are the relative 
frequency mismatch tA and the relative pulse amplitude (1 + 
/), see text. 
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FIG. 10: (color online) As in Fig. |pj but for the sym- 
metrized sequence BB[ CLJ) [(7r/2) 7r^(27r) 3 ^7r^(7r/2)o]. The 
pulse shapes are (a) lst-order pulses Si (the fidelity for the 
2nd-order pulses Qi is similar); (b) lst-order pulses with am- 
plitude correction; (c) 2nd-order pulses with amplitude cor- 
rection. Note how regular are the shapes of high-fidelity re- 
gions. 



five orders of magnitude with pulses Qi [Fig. He) . and 
by some two orders of magnitude for pulses Si [Fig. g(b)]. 
As for the Ising chain, the effect of the amplitude errors 
is weaker with specially-designed pulses L\ and L3 such 
that the lst-order coefficient v(f) scales as a higher power 
of /. With the pulse L\, the 2nd-order coefficient a is 
non-zero but small; it is seen from Fig. |l2|(b) that its ef- 
fect is to introduce a linear term 1 — F oc / which tends 
to skew the infidelity minimum away from / = 0. We 
should also note that with the Gaussian pulses G010 (not 
shown), even the on-resonance infidelity is out of range 
of the plots Fig. Ill 
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FIG. 11: (color online) As in Fig. |oj but for the sequence BUY 
[7ro7r,,i7r307r307r0] using only 7r-pulses. (a) lst-order pulses with 
amplitude correction; (b) 2nd-order pulses with amplitude 
correction. The absence of the error terms linear in a [see 
Eqs. (^), dSH)] produces a much wider high-fidelity region 
already with lst-order pulses. 

VI. CONCLUSIONS 

We presented a comprehensive study targeting pulse 
and sequence design and analysis based on a consistent 
high-order average Hamiltonian expansion. The numer- 
ical technique for expanding the evolution operator was 
originally introduced by us in Ref. and a complimen- 
tary analytical technique was developed for 7r-pulses by 
one of the authors in Ref. [2C]. 

The overall approach is to start with a closed sys- 
tem described by a finite-dimensional Hamiltonian Hs 
and design a sequence of shaped pulses such that the 
evolution operator would be accurate to a given or- 
der K in powers of Hs- The key to this approach 
are the NMR-style 1st- and 2nd-order self-refocusing 
one-dimensional pulses constructed for a single-qubit 
chemical-shift Hamiltonian (^). In this work we designed 
a number of such shapes for different rotation angles (f>o, 
and presented a careful analytical analysis of the first 
two leading orders of the average Hamiltonian theory 
for driven qubit evolution with the most general system 
Hamiltonian Hs- While any symmetric one-dimensional 
pulse shape is characterized by only three parameters, 
two of these can be set to zero by pulse shaping. The 
remaining parameter is also non-zero for an ideal "hard" 
((-function pulse. This leads to an important conclusion 
that the constructed pulses can be used as drop-in rc- 



FIG. 12: (color online) Decoupling errors as a function of 
relative pulse amplitude for a chain of n = 4 qubits with dif- 
ferent decoupling schemes as indicated, (a) Ising chain in the 
presence of individual chemical shifts Af. 8-pulse sequence 
with BBi composite pulses offers the best accuracy, which re- 
mains essentially the same whether the sequence is used with 
regular 1st- or 2nd- order pulses or with the pulses stabilized 
against amplitude errors (lst-order pulses L\ and 2nd-ordcr 
pulses 1/3 ). However, the pulse L3 works well enough even 
with regular 8-pulse sequence. While the details of the ampli- 
tude scaling differ, at the level oil — F — 10 -4 , the lst-order 
amplitude-protected pulses L\ and regular 2nd-order pulses 
Q\ have comparable accuracy. The use of lst-order pulses 
show relatively poor performance even on resonance, (b) XXZ 
chain in the presence of individual chemical shifts Af . With 
XXZ coupling, The BBi composite pulse is no longer accu- 
rate, as the errors appear already in the linear order in r J x 
(not shown). With the regular 8-pulse sequence 8, the best 
accuracy is obtained for the pulses with amplitude correction. 



placement for hard pulses; with proper pulse placement 
the results should be identical to first two orders. The 
structure of errors appearing in higher orders of the evo- 
lution operator can be understood by analyzing the nu- 
merical time-dependent perturbation series for the evo- 
lution operator of a closed system. 

An important advantage of this approach is that the 
expansion order offers a natural classification of the er- 
ror operators. As a result, (i) the convergence regions 
have regular shapes as a function of parameters [see 
Figs. |lO](c,d) and |ll](b,c)]. Furthermore, with local two- 
(or few-)qubit couplings dominant, (ii) the error opera- 
tors can be placed on connected clusters of up to k + 1 
qubits for terms of order k 7 which allows one to under- 
stand their structure in terms, e.g., the direct products of 
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up to (k + 1) Pauli matrices. Once their structure under- 
stood, the convergence can be readily improved by sup- 
pressing the error operators, as in our analysis of pulse- 
amplitude errors. We emphasize, that such an analy- 
sis can be performed even for very large qubit systems. 
Thus, (iii) this approach is characterized by scalability 
with the system size, as we illustrated by analyzing de- 
coupling infidelity with the system size [Fig. § . Although 
in this work we concentrated on the dynamics of closed 
systems, another important advantage is that (iv) the 
high-order control sequences result in lower deceherence 
in the presence of slow environmental modesB3Ej. 

Most obvious application of highly-optimized shaped 
pulses of the sort presented in this work is in solid-state 
quantum computation, where the bandwidth available 
for quantum gates is typically limited. Our techniques 
based on analytical and numerical high-order average 
Hamiltonian theory offers a systematic scalable approach 
for constructing gates for such multi-qubit systems, with- 
out need of solving their full dynamics. However, even if 
the bandwidth does not appear to be at premium, simple 
pulse shaping (e.g., using lst-order pulses) can still offer 
a substantial improvement of control accuracy. 
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APPENDIX A: AVERAGE FIDELITY 

Here we discuss the calculation of the fidelity averaged 
over the initial state, in the case of unitary evolution with 



known evolution matrix U, while the desired evolution 
matrix is Uq. Let us write the density matrix of the initial 
pure state as po = iptp where ip is an iV-component 
complex vector. Then the actual density matrix is p = 
UijjipW^ , while the desired density matrix is pidcai = 
U ^U&. The fidelity with the given initial state 



F,i, 



tr(p ideal p) = tv{U ^ulu^U^) 



(Al) 



ijkl 



The only condition on the components ipi of the wave- 
function is the normalization, 1 = ip^tp = J^. \ipi\ 2 - Gen- 
erally, this means that the average of the product in 
Eq. QAl| ) can only depend on the identity tensor Stj . By 
symmetry, (ipiip^ip k ipf) = A{8 lj 5 k i + 8 a 8 k j), where the 
unknown coefficient A can be computed from the nor- 
malization by tracing over i = j, k = I. We obtain 
1 = A(N 2 + N), so that the average fidelity 



F 



N+ \trV\ 2 
N + N 2 



l-F = 



N 



trtf 



N + N 2 



(A2) 



where V = UqU. Numerically, with V close to identity 
matrix, the loss of precision can be avoided by expressing 
the infidelity 1 — F in terms of the modified mismatch, 



<5 2 = tr(l-F) t (l-F), V = V 



tvV\ 
try 



(A3) 



Namely, since 5 2 = 2N — 2|trV|, the average infidelity 
can be written as 



1 - F 



S 2 (AN-S 2 ) 
A(N + N 2 )' 



(A4) 
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